function [Z, Zprob] = tauchen(N, mu, rho, sigma, m)
    Z = zeros(N,1);
    Z(end) = m * sqrt(sigma^2 / (1-rho^2));
    Z(1) = -Z(end);
    
    zstep = (Z(end) - Z(1)) / (N-1);
    
    Zprob = zeros(N, N);
    a = (1-rho)*mu;
    
    Z(2:end-1) = Z(1) + zstep * (1:N-2);
    Z = Z + a / (1 - rho);
    
    for k = 1:N
        if k == 1
            Zprob(:,k) = normcdf( (Z(1) - a - rho * Z + zstep / 2) / sigma ) ;
        elseif k == N
            Zprob(:,k) = 1 - normcdf((Z(N) - a - rho * Z - zstep / 2) / sigma);
        else
            Zprob(:,k) = normcdf((Z(k) - a - rho * Z + zstep / 2) / sigma) -...
                      normcdf((Z(k) - a - rho * Z - zstep / 2) / sigma);
        end
    end
    
end
    